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ABSTRACT 

We present a stellar dynamical estimate of the black hole (BH) mass in the Seyfert 1 
galaxy, NGC 4151. We analyze ground-based spectroscopy as well as imaging data 
from the ground and space, and we construct 3-integral axisymmetric models in order 
to constrain the BH mass and mass-to-light ratio. The dynamical models depend on 
the assumed inclination of the kinematic symmetry axis of the stellar bulge. In the 
case where the bulge is assumed to be viewed edge-on, the kinematical data give only 
an upper limit to the mass of the BH of ~ 4 x 10^ Mq (Ict). If the bulge kinematic 
axis is assumed to have the same inclination as the symmetry axis of the large-scale 
galaxy disk (i.e., 23° relative to the line of sight), a best-fit dynamical mass between 
4 — 5 X 10^ Mq is obtained. However, because of the poor quality of the fit when the 
bulge is assumed to be inclined (as determined by the noisiness of the surface and 
its minimum value), and because we lack spectroscopic data that clearly resolves the 
BH sphere of infiuence, we consider our measurements to be tentative estimates of the 
dynamical BH mass. With this preliminary result, NGC 4151 is now among the small 
sample of galaxies in which the BH mass has been constrained from two independent 
techniques, and the mass values we find for both bulge inclinations are in reasonable 
agreement with the recent estimate from reverberation mapping (4.57l'^Q4y ^ ^0) 
published by Bentz et al. 

Subject headings: black hole physics — galaxies: active — galaxies: nuclei — galaxies: structure 
— stellar dynamics 
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1. Introduction 

Evidence for the existence of supermas- 
sive black holes (BHs) in galactic centers 
has progressively become stronger over the 
past decade (see iFerrarese fc FordI l2005l . for 
a recent review). Direct kinematical mea- 
surements of BH masses (Mbh) from trac- 
ers of various types have been made in 
roughly three dozen galaxies. This growth 
in the number of BHs with mass measure- 
ments has led to the discovery of corre- 
lations between Mbh and various proper- 
ties of the host galaxies. The measured 
masses have been fou nd to be correlated with 
bulge luminositv (IKormendv &: Richstone 
1995:'M cLure k Dunlop ll2002l : lMarconi fc Hunt 



(AGNs). In such systems, the brightness 
of the non-stellar continuum flux from the 
AGN complicates direct measurement of the 
BH mass from kinematical tracers around 
the central object. However, BH masses 
measured in broa d-line AGNs using rever- 
beration mapp ing (jBlandford McKedfToii : 
PetersonI Eooi ) have also been found to foL 



low a Mrh — rela ti on (iGebhardt et al 



2000bl : IFerrarese et al.1 l200ll : lOnken et al 



2OO3I : iNelson et al.l |2004| V Reverberation 



mapping relates the properties of gas in 
the broad-line region (BLR) to variability in 
the ionizing continuum flux from the central 
AGN. The time delay between variations in 
the continuum and the response of a broad 

II ! ■ " — -■ — I ' eisission line arising from the BLR is used to 

2003|;iErwin et al."2004'), bulge massjMerrit t^fc Ferrag^ ^mine the characteristic radius (r) of the 

BLR, and the width of the emission line is 
used to estimate the line-of-sight velocity dis- 



RiJ 

2004 ) . the central velocity dispersion of the 



2001: iMarconi fc Hunt 200^: |Harme 



Ferrarese &: Merritt 


2OO0F 


Gebhardt et al. 


2000a), the degree of concentration of bulge 


light ( 


Graham et al. 


2OO1I: Erwin et al. 


2004; 



Graham &: Driverl 120071 ). and, possibly, the 



mass of the su rrounding dark matter halo 
(|Ferraresel 12003 ). 

Most of the BHs that were used to de- 
fine the Mbh — c* and other scaling relation- 
ships between BH masses and host galaxy 
properties reside in quiescent galaxies. How- 
ever, BHs also power active galactic nuclei 
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persio n of the gas at that radius. I Collin et al 
(I2OO6I ) showed that measuring the line width 
as the second moment of the variable part 
of the emission line profile (aune) is less bi- 
ased than using the full width at half max- 
imum (FWHM). The mass is then given as 
Mbh = /^crfine/^' where G is the gravita- 
tional constant and / is a scaling factor that 
depends on the BLR geometry and kinemat- 
ics. Because the BLR is not spatially resolved 
with current technology, mass measurements 
from reverberation mapping (and, specifically, 
the value of /) must be calibrated against 
BH masses derived through other methods. 
One way to determine / is to use the cor- 
relation between central stellar velocity dis- 
persion (fj*) and A/bh for those galaxies in 
which rev erberation ni a ,pping estimates are 
available. lOnken et al.l ( 20041 ) calibrated the 
Mbh — o"* relationship for 16 reverberation- 
mapped AGNs and calculated the average 
scaling factor (/) = 5.5 it 1.8. 

Another method of calibrating the BH 
masses from reverberation mapping is to ap- 
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ply one of the techniques commonly used for 
estimating BH masses in quiescent galaxies 
(e.g., dynamical mass estimation from stellar 
kinematical data) to a reverberation-mapped 
AGN. This approach is an important step in 
checking the masses of BHs determined from 
reverberation mapping since it does not rely 
on a secondary indicator of the mass (like cr* 
in the case above). To firmly constrain BH 
masses, the kinematic measurements must 
reach inside the sphere of influence, the re- 
gion wherein the stellar dynamics are domi- 
nated by the BH. Satisfying this requirement 
is particularly challenging in AGNs, in which, 
at the small angular separations correspond- 
ing to the radius of the sphere of influence 
(rh = GMbh/c^), the stellar light is eas- 
ily overwhelmed by the luminosity from the 
AGN. 

Based on the revised BH masses pro- 
vided by the reve rberation reanalysis of 
Peterson et al.l (j2004l ) and correspo nding bulge 



veloc i ty dispersions measurements (|F errares e et al 



2OO1I : lOnken et all l2004l : iNelson et al. J004), 
there are two reverberation-mapped AGNs 
for which the BH sphere of influence should 
be spatially resolvable with existing facili- 
ties: NGC 3227 (as mentioned above) and 
NGC 4151, both local broad-line Seyfert 
galaxies. The BH masses in these two AGNs 
were recently estimated through observations 
of their nuclear gas kinematics and by mod- 
eling the gas velocity field as a flat, circu- 
lar disk dHicks k Malka^] l2007l). In addition, 
Davies et al.l (j2006l ) measured the BH mass 
in NGC 3227 from stellar dynamical model- 
ing. We discuss these results further in § [5l 
A handful of narrow-line AGNs also have dy- 
namical estimates of their BH masses, but 
because narrow-line objects are not capable 
of being reverberation mapped, these objects 
cannot act as calibrators between the two 
mass scales. 

Here we present the results of a program to 



measure the BH mass in NGC 4151 via stel- 
lar kinematical modeling. In § [2l we describe 
the imaging and spectroscopic observations as 
well as the data analysis process. In § [3l we 
provide a brief overview of the standard 3- 
integral orbital superposition method of stel- 
lar dynamical modeling that we use to fit the 
data. In § m we present the results of the 
modeling. In § O we discuss the implications 
of our results for the calibration of BH masses 
from reverberation mapping and future steps. 

2. Observations and data analysis 

2.1. NGC 4151 

NGC 4151 is a disk galaxy with faint spiral 
arms, a weak large-scale bar, and a promi- 
nent bulge. It is one o f the or i ginal s ample of 
galaxies identified by Sevfert ( 19431 ) as hav- 
ing unusual nuclear emission line properties, 
and has been studied exte nsively over a wide 
range of wavelengths (see Ulrich 20001 ). The 
heliocentric r adial velocity of t he galaxy is 



998 km s (jPedlar et al.l Il992l ) , implying a 



distance of 13.9 Mpc if in pure Hubble expan- 
sion with Hq = 72 km Mpc~^. However, 
other estimates of the distance to NGC 4151 
have ranged from 10 to 30 Mpc based on var- 
ious estimates of the Virgoce ntric infall cor- 



rectio n (see the discussion by iMundell et al 
I999I ). We adopted a distance of 13.9 Mpc, 



but note that our dynamical estimate of the 
BH mass scales linearly with distance. For 
ex ample, apply i ng th e local infall correction 
of lMould et aD (|200d ) would give a distance 
of 17 Mpc, increasing our estimated masses 
by roughly 20%. At our assumed distance of 
13.9 Mpc, 1" corresponds to 67 pc. 

The (distance-independent) reve rberation 

mass estimate for NGC 4151 given bv lPeterson et al 
(|2004l ) was 1.33x10^ M© (from a combina- 
tion of Ha and H/3 measurements). How- 
ever, a recent reanalysis of archival UV mon- 
itoring data from the International Ultravi- 
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olet Explorer (jMetzroth et alj 1200611 and a 



new H /3 reverberation campaign (|Bentz et al 



2006bl l indicate a higher BH mass: 4.57+|];|^ 



10' Mq. NGC 4151 has a bulge velocity 
dispe rsion of 93ibl4 km s~^ (jFerrarese et al 



200 ll ). which, when combined with the Bentz 
et al. reverberation mass, gives a sphere of 
influence with rh = 23 pc or 0'.'34. 

2.2. Imaging Data 

We made use of a set of ground-based B VR 
images from the MDM 1.3 m McGraw-Hill 
telescope as well as Hubble Space Telescope 
(HST) imaging data obtained with the High 
Resolution Channel (HRC) of the Advanced 
Camera for Surveys (ACS). Both datasets 
are part of an independent project designed 
to study host galaxy conta mination in AGN 
spectra (jBentz et al.ll2006al ). 



2.2.1. Imaging with HST 

The ACS/HRC data were obtained with 
the F550M filter and have a total exposure 
time of 1020 s. The full details of the data re - 
duction are described by Bentz et al. ( 2006al ). 
and the final HST image is shown in Figured) 
The image has a plate scale of 0'.'025 pixel"^ 
and a field of view of ~ 25" x 29", with a point 
spread function (PSF) FWHM of ^'06. 

2.2.2. Ground-based imaging 

The MDM 1.3 m McGraw-Hih telescope 
was used with the "Templeton" CCD to take 
direct images of NGC 4151 in BVR. The re- 
spective total exposure times were 1470, 1200, 
and 1220 s. The images in each filter were 
bias-subtracted, flat-fielded, registered, and 
then median-combined. The i?-band image 
(Fig. [2]) was used to determine the surface 
brightness proflle on large scales (^ 12. 3p . The 
image has a plate scale of 0'.'508 pixel"^, a 
field of view of roughly 8' on a side, and a see- 
ing FWHM of « 2"3. Flux calibration was 



achieved through aperture photometry of a 
star in the image ( Eyermann et al. 20051 ). 
Beh fc de .TomI (|200lh employed models oi' 



stellar population synthesis and spiral galaxy 
evolution, in conjunction with measured rota- 
tion curves, to derive relationships between 
galaxy color and stellar mass-to-light ratio. 
We used our multi-color imaging data to es- 
timate the i?-band stellar mass-to-light ratio, 
Tr (in units of solar masses per i?-band solar 
luminosity throughout), as follows. The im- 
ages were aligned, convolved with Gaussians 
of appropriate size to bring the PSFs to the 
same size as the band with the worst seeing, 
and a pixel-by-pixel color map was created 
(Fig.E]). The (B-R) and (V-R) colors were 
then determined ctS Si function of radius. The 
strong decrease in {B — R) and (V — R) near 
the center is consistent with AGN contami- 
nation, therefore we ignored the central 3". 
With (5-i?)wl.5 mag and ( V- R)^0.6 mag 
we e s timated Tp to be « 3 (IBell fc: de Jong 
200ll : iBeU et al.llioO^ ) . 



2.3. Surface Brightness Profile Model- 
ing 

We derived the surface brightness profile 
from the MDM i2-band image and the HST 
F550M image. As the goal of this procedure 
was to develop a model for the mass distribu- 
tion in the galaxy, we needed to remove the 
AGN contribution to the light profile and ex- 
tract the stellar surface brightness. The pro- 
cess of removing the AGN light was compli- 
cated by the very different spatial character- 
istics and filter bandpasses of the two images. 
We used GALFIT0, the p ublicly available 2 - 



D galaxy fitting software (jPeng et al.l l2002l ) 



to remove the AGN contamination. In order 
to most accurately subtract the AGN fiux, we 
performed a multi-component fit to the major 
surface brightness features of the galaxy and 



^http: //zwicky . as . ajrizona. edu/$\sini$cyp/work/galf it/galf it .html 
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jointly determined the galaxy properties and 
the AGN flux to be removed. 

For the MDM R-hand image, we sub- 
tracted a sky level equal to the median flux 
beyond the visible extent of NGC 4151, and 
fit two Sersic profiles (one for the disk and one 
for the bulge), allowing the Sersic index to be 
a free parameter in each component. We also 
used a PSF for the AGN emission that was 
based on one of the stellar profiles in the im- 
age. The fit was done in several steps, initially 
masking out the bulge and nucleus in order 
to fit the disk. Next, we removed the mask 
from the bulge and used the earlier results 
from the disk-only fit as the starting values 
for the disk component. Then we removed 
the mask from the nucleus and allowed the 
nuclear magnitude to vary. We left the mask 
in place over several stars in the field and 
over a number of bright knots near the ends 
of the weak bar in NGC 4151. Because of 
slight differences in the profile shapes of the 
stehar PSF and the AGN light (possibly due 
to small image distortions between the posi- 
tions of the AGN and PSF star in the field of 
view), manual adjustment of the PSF mag- 
nitude and position was required to produce 
the final fit. The best-fit Sersic index values 
for the disk and bulge components were 0.73 
and 1.48, respectively. 

The small field of view of the HST im- 
age meant there were no galaxy-free regions 
from which to measure the sky fiux. In- 
stead, the background count rate for the HRC 
in the F550M filter was retrieved from the 



ACS Instrument Handbook ( Pavlovsky et al 



20041 ) and was multiplied by the total expo- 
sure time to estimate the background flux. 
The effective radii of the two Sersic profiles 
from the fit to the MDM image were con- 
verted to the appropriate pixel scale for the 
HST image and held fixed. The total mag- 
nitudes of the components were allowed to 
vary to account for differences arising from 



the different bandpasses between the images. 
We attempted to model the HRC PSF with 
Tiny Tir4j, which produces PSFs for the var- 
ious HST instruments at any desired posi- 
tion in the field. Tiny Tim also allows one 
to input a user-defined spectrum to better 
match any chromatic distortions in the op- 
tical path. We joined two archival Space 
Telescope Imaging Spectrograph (STIS) spec- 
tra of NGC 4151 taken with the G430L and 
G75 0L low-resolution g ratings and a 52" x O'.'l 
slit (jNelson et al.ll200d ^. and fed the resulting 
spectrum to Tiny Tim. However, the output 
PSF proved to be a poor match to the nuclear 
emission. We searched for all archival stel- 
lar images taken with the HRC in the F550M 
filter. Unfortunately, this search turned up 
only a handful of white dwarf observations, 
all with relatively short exposure times. The 
best of these was a 2 s im age of BD+17°4708 
(|Bohhn fc Gilhlandlliooi ). and this image was 
then used as our HRC PSF. 

There were some strong residuals around 
the best GALFIT fit to the HST image, but 
it was unclear whether they were primarily 
due to dust obscuration or to [O III] A5007 
emission leaking into the F550M bandpass. 
To assess the latter possibility, we retrieved 
an [O III] image of NGC 4151 taken with the 
linear ramp filter on the Wide Field Plane- 
tary Camera 2 (WFPC 2) instrument of HST 
(|Hutchings et al.lll999l ). After registering the 
images and resampling the HRC image to the 
slightly larger pixel scale of WFPC2, we were 
unable to find a scaling factor for the [O III] 
image that produced a satisfactory subtrac- 
tion of the remaining residuals. While some 
of the positive-flux regions in the residual im- 
age may have been due to [O III] emission, it 
is likely that dust obscuration also played an 
important role in shaping the F550M image 
morphology. Ultimately, these regions were 



^ ihttp: //www, stsci . edu/sof tware/tinytim/tinytini.htmll 
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simply included as-is in the GALFIT fitting 
procedure. 

With the AGN emission subtracted, we em- 
ployed the Multi-Ga ussian Expansion (MGE) 
routine described by ICappellaril (|2002l ) to de- 
project the surface brightness distribution 
into a 3-D luminosity model. The MGE rou- 
tine fits a series of ~ 10 concentric elliptical 
Gaussians to the 2-D galaxy image. (The 
program has been made publicly availabl^ 
for use with IDlS.) The fit was performed si- 
multaneously to the MDM and HST images, 
and this allowed us to rescale the HST flux 
level to match that of the MDM image where 
the two overlapped in radius. The scaling 
was consistent with that expected from the 
differences in plate scale and filter bandpass. 
We then tied the overall flux calibration to 
that determined by the stellar photometry 
in the i?-band. This approach had the ad- 
vantage that the dynamical Tr, values could 
then be compared to the Tr estimates from 
the multi-color MDM photometry. Adopt- 
ing an i?-band extinctior0 in t he direction of 



NGC 4151 of ^i?=0.074 mag (jSchlegel et al 



19981 ). we converted the MGE output of flux 
in each Gaussian to solar luminosities (Ta- 
bled]). The 2-D contours of the MGE model 
are shown in Figure U] (with the same size 
scale as Figures [U & [2|), and the major-axis 
i?-band surface brightness profiles of both the 
raw data (determined with the IRAF task EL- 
LIPSE) and the AGN-free MGE output are 
shown in Figure O 

2.4. Spectroscopic Data 

We acquired spectroscopic data at the 
MMT Observatory, and we revisited the Kitt 
Peak National Observatory (KPNO) spec- 



http: //www. strw. leidenuniv.nl/$\sim$mcappell/idl/ 



' jhttp: //www.rsinc . com/idl/| 
® Retrieved from 
|http : / /nedwww . ipac . caltech . edu/ 1 



tra of NGC 4 1 51 th at had been used by 
Ferrarese et al. (|200lh to determine a*. The 
stellar dynamics were determined using the 
Ca II triplet (hereafter CaT) stellar absorp- 
tion lines at AA8498, 8542, and 8662 A. We 
extracted the stellar kinematics along the two 
(different) position angles (PAs) of the MMT 
and KPNO data. 

2.4.1. MMT Spectra 

We observed NGC 4151 and velocity stan- 
dard stars with the Blue Channel Spectro- 
graph at the 6.5 m MMT on UT 29 May 
2004. We used the 1200 lines mm^^ grating 
with the LP-530 order-blocking filter, a l" slit 
width, and the "ccd35" detector. The spectra 
were centered at 8600 A for measurement of 
the CaT with a dispersion of 0.50 A pixel"^ 
and a resolution of 1.4 A. The plate scale 
was 0'.'3 pixe^^ and the PA was 69°. While 
the photometric conditions were quite good, 
strong and gusty winds caused both poor see- 
ing and problematic guiding (yielding a final 
spatial profile with FWHM~ 3"). Our total 
exposure time was 11700 s. 

The spectra were rectified and wavelength- 
calibrated using the bright sky lines that per- 
meate the observed spectral region (the same 
wavelength solutions were applied to the spec- 
tra of a K-giant star, HR 4521, as the stellar 
exposures were too short for sky lines to ap- 
pear and the lines present in the HeNeAr cali- 
bration lamp spectra were much more sparse) . 
Flat field spectra taken among the obser- 
vations of NGC 4151 were employed to re- 
move the strong fringing (20% peak-to-peak) 
present in the spectra. Although the correc- 
tions left residual fringes of nearly 10%, the 
spectral regions around the bluest two lines 
of the CaT were free from contamination, and 
so we limited our subsequent analysis to those 
lines. 



NED: 
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24-2. KPNO Spectra 



The work of iFerrarese et alj (|200lh used 
the Ritchey-Chretien Spectrograph on the 
Mayah 4 m telescope at KPNO on UT 9 April 
2001 to measure the bulge velocity dispersion 
ofNGC4151. Det ails of the observations were 
fully described in Ferrarese et al. ( 200ll ). but 
we restate the key features here. The obser- 
vations were taken with a 2" slit at PA=135°, 
the detector plate scale was 0'.'69 pixel"^, and 
the spectral resolution was ~ 1.7 A. The total 
exposure time was 3600 s. 

The observations were rectified and wavelength- 
calibrated with the sky lines present in the 
spectra, and then combined. Spectra of the 
star HR 4521 were obtained for the absorption 
line profile measurements and were reduced in 
the same way as the AGN data. 

2.5. Spectroscopic Analysis 

The MMT spectra were extracted in 1-pixel 
(O'.'S) apertures out to ±12" (±40 pixels), 
tracing along the peak of the spatial profile. 
While the seeing disk (modeled Gaussian 
with fj = 1'.'27) was much larger than our in- 
dividual extraction windows, the dynamical 
modeling routine took the seeing into account 
when fitting to the observed kinematics. To 
isolate the CaT lines for profile fitting, the ex- 
tracted spectra were normalized with a high- 
order spline fit to the continuum, any remain- 
ing fringing, and several AGN emission fea- 
tures. 

Because of the lower signal-to-noise ratio, 
S/N, in the KPNO data, extractions were 
made with a variable window size, ranging 
from 1 pixel (0'.'69) at the peak of the light to 
a 7-pixel width in the furthest bins (extend- 
ing to distances of ~20 pixels, or 14"). The 
seeing for these observations was modeled as 
a Gaussian with a = 0'.'76. 

The spectroscopic data were analyzed with 
the Penalized Pixel-Fitting (pPXF) method 



of ICappellaril (|2004l rl. The pPXF program 
finds the best fit to the line-of-sight veloc- 
ity distribution (LOSVD) based on a max- 
imum penalized likelihood approach to the 
broadening of a template spectrum. The 
method parametrizes the LOSVD in the 
stan dard way as a Gauss - Hermite serie s 
der Marel fc FranxlEggj : iGerhardlEggj ) , 



van 



which we limited in our analysis to the first 
four terms {V, a, h^, h^), where V is the ra- 
dial velocity, a is the velocity dispersion, and 
hs and /14 are the coefficients of the next two 
elements of the Gauss-Hermite expansion. 

The pPXF routine uses a stellar spectrum 
as its template and simultaneously fits to 
all desired Gauss-Hermite terms over a user- 
defined wavelength region. The ability to 
select multiple fitting regions allowed us to 
exclude regions of noise or AGN emission 
that might lie between the stellar absorp- 
tion lines we are trying to measure. Each 
of the extracted spectra from the MMT and 
KPNO data were run through the pPXF algo- 
rithm with a corresponding stellar template. 
HR 4521 was the only template common to 
both datasets. Thus, to minimize systematic 
differences between the MMT and KPNO re- 
sults, and because the CaT is generally in- 
sensitive to pr oblems froni mism atched stel- 
lar templates ( Barth et al. 20021 ). we limited 
our analysis to the single template star. The 
resulting kinematic parameters are shown in 
Figure M and listed in Table [21 

In addition to the narrow-aperture extrac- 
tions described above, the MMT spectrum 
was extracted with a width of ~ 4" to provide 
an independent measurement of the bulge ve- 
locity dispersion. Applying the pPXF proce- 
dure gave cj* = 109 ± 11 km s~^, reasonably 
consistent with the results of both Ferrarese 
et al. (2001; 93 ± 14 km'^) and Nelson et al. 



^The IDL routine is publicly available from 
http: //www, strw. leideiiuiiiv.iil/$\sim$mcappell/ idl/[ 
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(2004; 97 ±3 km s"^). 

2.6. Target-of-Opportunity Observa- 
tions with STIS 

One way to minimize the problem of nu- 
clear glare when attempting to measure stel- 
lar features in an AGN is to take advantage 
of the variable nature of the nuclear emis- 
sion and observe the target when the AGN 
is dim. We obtained Target-of-Opportunity 
(ToO) observations of NGC 4151 with STIS 
on HST. Monitoring of both NGC 4151 and 
NGC 3227 was carried out with ground-based 
spectroscopy at the Crimean Astrophysical 
Observatory (CrAO) and with space-based X- 
ray flux measurements from the Rossi X-ray 
Timing Explorer (RXTE). 

The threshold for triggering the ToO obser- 
vation had been derived from ground-based 
i?-band observations and bulge-disk-nucleus 
decompositions of NGC 3227 and NGC 4151 
( Virani et al. 2000l ) . These data allowed us to 
estimate the stellar flux we could expect to fall 
within the STIS slit and also to independently 
examine the effects on the S/N of varying the 
nuclear brightness. We assumed a stellar con- 
tribution near the center of each galaxy to be 
roughly like an elliptical galaxy (flat in Fx), 
with a spatial distributi on given by the bulge 
fit of Virani et al. (|2000l ). Thus, we were able 
to parametrize the exposure time as a func- 
tion of the nuclear non-stellar flux. We cal- 
culated the thresholds of nuclear flux for each 
target that allowed us to reach our desired 
stellar S/N of 50 A^^ in a reasonable amount 
of exposure time. For NGC 4151, this thresh- 
old was equal to 4 x 10~^^ erg cm~^ A.~^, 
which the AGN had been near in the recent 
past. 

On 2003 November 19, the CrAO obser- 
vations of NGC 4151 indicated that the con- 
tinuum flux level had fallen below the trig- 
gering threshold. Figure [7] shows a historical 
light curve for NGC 4151 and indicates with 



arrows the date of a previous STIS observa- 
tion (a shorter exposure that employed the 
STIS occulting bar to block the nuclear light, 
but which also blocked the stellar signature 
near the BH), the date when our ToO obser- 
vation began, and when the STIS instrument 
stopped functioning. 

The STIS observations were executed on 
UT 2003 December 15 and 19. The instru- 
mental setup consisted of the G750M grating 
and the 52" x O'.'l slit, oriented with a PA of 
69°. The total exposure time was 29730 s. 
In addition to the observations of NGC 4151, 
spectra of HR 4521 were obtained on UT 
2003 December 26 with the same instrumen- 
tal setup to act as our comparison for the CaT 
line profiles. 

While the continuum flux stayed close to 
the threshold value between the ToO trigger- 
ing on 2003 November 20 and the actual time 
of observation, the extracted STIS spectrum 
does not show the presence of the CaT lines, 
and so we were unable to use the STIS data to 
constrain the stellar kinematics near the BH 
in NGC 4151. 

The lack of absorption lines can be used 
to place upper limits on the number of gi- 
ant stars in the central region of NGC 4151. 
The spectrum of the comparison star HR 4521 
was used to estimate the maximum number 
of such stars that could exist in the center 
of NGC 4151 without producing measurable 
lines. HR 4521 has a Hipparcos parallax mea- 
surement of 15.80ib0.59 mas, implying a dis- 
tance of 63.3ib2.4 pc. We then scaled the 
luminosity of HR 4521 to the flux appropri- 
ate for our assumed distance of NGC 4151 
(13.9 Mpc), and found the multiplicative fac- 
tors required to produce visible spectral fea- 
tures when subtracting the scaled stellar spec- 
tra from NGC 4151 spectra of different aper- 
ture sizes. For the inner O'.'l (6.7 pc), we de- 
rived an upper limit of 7.5 x 10^ giant stars. 
A wide-aperture extraction (~ 3" x O'.'l, or 
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~200 pc X 6.7 pc) only increased the up- 
per limit to 1.25x10^ stars. For comparison, 
the mass enclosed within a radius of 7 pc at 
the center of the Milky Way is ~ 2 x 10^ 
Mq, only about 10% of which is acco unted 
for by the BH (see lMeha FalckelboOlK Ap- 
proximately ~1% of the stellar mass within 
that region of t he Milky Way is c omposed 



is the one described by Valluri et al. (2004; 
hereafter VME04) and we refer the reader to 
that paper for a detailed description of the 
method. 

We constructed a 3-D mass model of the 
galaxy by deprojecting the surface bright- 
ness distribution produced by the MGE rou- 
tine 12.31), usin g the method outlined by 



of bright giants |Genzel et alJ|l996|), giving ICappellaril (l2002l ). and assuming an error of 



roughly 2x10^ M© and the equivalent number 
of stars (assuming an old bulge stellar popula- 
tion). Whether we take the NGC 4151 bulge 
mass estimate of Wandel (2002; a factor of ~4 
smaller than the Milky Way bulge mass given 
by Bissantz et al. 1997), or we estimate the 
bulge mass from the MBH-A^buige relation of 
Haring & Rix (2004; a factor of 3 larger than 
the Milky Way bulge), the lack of CaT ab- 
sorption is consistent with NGC 4151 having 
a similar bright giant fraction as the Milky 
Way. 

As noted bv iBentz et"al] (|2006al ). ground- 
based imaging data often are insufficient 
to make an accurate galaxy decomposition. 
Thus, our reli ance on the g round-based imag- 
ing analysis of lVirani et al.l (H) to estimate 
the appropriate exposure time for the STIS 
observation may have contributed to our fail- 
ure to detect the stellar absorption features. 
However, it is also possible that the CaT lines 
were filled in by emission arising from either 
young stars near the center of NGC 4151 or 
from the AGN itself. 

3. Dynamical Modeling 

We carried out dynamical modeling of the 
kinematical and imaging data using the now 
standard met hod of orbit s u perpo sition first 
described by Schwarzschild ( 19791 ). This al- 
lowed us to set dynamical constraints on the 
possible values of the BH mass (Mbh) and the 
stellar mass-to-light ratio (Tr) that are con- 
sistent with the surface brightness and spec- 
tral data. The modeling algorithm we used 



10~^ in each mass cell. We also assumed an 
inclination angle of the bulge (described be- 
low) and, to normalize the mass density, a 
value for Tr, which was constant with radius. 
An unknown point mass (Mbh) representing 
the central supermassive BH was added to 
the mass distribution and a large set of orbits 
were integrated in the resultant gravitational 
potential for every pair of parameters (Mbh; 
Tr), building an orbit library for each (Mbhj 
Tr) combination. We constructed models for 
a wide range of values of Mbh and Tr (12- 
17 different values of Mbh and 8-16 values of 
Tr). 

Two different inclination angles were used 
to model the bulge. The large-scale galaxy 
disk is somewhat inclined to the line of sight 
(the rotation axis of the disk having an in- 
clination of_J_^^ 20_— 25° relative to the line 
of sight; ISimkinI Il975l ). However, the bulge 



appears very close to circular in projection 
(axial ratio q ~ 0.99), suggesting that it is 
nearly spherical. We therefore constructed 
two sets of models. We first assumed that the 
bulge is viewed edge-on (i = 90°) and has its 
kinematic symmetry axis in the plane of the 
sky. The assumption of edge-on axisymmetry 
should be reasonable for orbits near the BH. 
We also constructed a second set of models 
for which we assumed that the rotation axis 
of the axisymmetric bulge has an inclination 
of i = 23°. Because of these motivations for 
the two particular values of the inclination, 
we treat the cases independently, rather than 
marginalizing over the inclination. 
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The presence of a weak, large-scale bar in 
NGC 4151 is noticeable beyond 50" and is 
expected to result in non-axisymmetric mo- 
tions in the stellar orbits, which cannot be 
modeled by our axisymmetric code. Although 
the bar is expected to b e kinematically weak 
(jMundeh fc Shond Il999l ). we restricted our 
models to fit the light distribution between 
O'/OOOIS (0.01 pc) and 50" (3.35 kpc). Again, 
within this radius, the projected isophotes 
were almost circular. 

The total number of luminosity constraints 
from the surface brightness profile was 266 
(the inner 196 from the HST photometry, and 
the remainder from the MDM photometry). 
Kinematical constraints (V, o", /i3, /14) were 
available for 58 apertures within ~ 15" from 
the KPNO and MMT observations. Since the 
kinematical quantities (especially V and a) 
obtained with both spectrographs were not 
symmetrical with respect to the center of the 
galaxy, the data at each slit position +x were 
averaged with those from X , clS IS standard 
procedure. The Schwarzschild orbit superpo- 
sition code assumes that the mass distribution 
is axisymmetrical. Consequently when the 
mass contribution of an orbit is stored in the 
library, it is stored on a grid in the meridional 
plane of the galaxy in cylindrical polar coor- 
dinates. However since we fit the kinematical 
data in the plane of the sky, the kinematical 
information from the orbit library was stored 
in the plane of the sky. It is necessary to sym- 
metrize the kinematical data because if one 
does not do so, the model will try to simul- 
taneously fit non-axisymmetric kinematical 
features and an axisymmetric mass distribu- 
tion. Since these conditions are physically 
inconsistent with each other, the solutions 
obtained have poor and are biased. The 
total number of kinematical constraints mod- 
eled were 4 Gauss-Hermite parameters x 29 
apertures= 116. 

Each model was constructed with an or- 



bit library consisting of 8100 different orbits. 
The contributions of each orbit to the ob- 
served LOSVDs and the deprojected lumi- 
nosity distribution were determined. A non- 
negative least squares p rogramming routine 
(|Lawson &: Hansonlll995l ) was used to find the 
weighted superposition of the orbits that best 
reproduced both the assumed model stellar 
density distribution and the observed kine- 
matical data. 

Smoothing, or "regularization" , of the or- 
bital solutions is typically performed in this 
type of analysis, and a variety of methods 
have been employed to that en d. One can lo- 



cally smooth the phase space (jCretton et al, 
I999I ). i mpose a "maximum entropy" con- 



straint dCebhardt et al.1 l2003l ^. or include a 
penalty function with an adjustable smooth- 
ing parameter (VME04). In the case of the 
penalty function, poorly selected values of 
the smoothing parameter can generate mis- 
leading results (VME04). For example, over- 
smoothing effectively restricts the range of al- 
lowed orbits, and can give erroneous best-fit 
values of Mbh- Here we present models with- 
out regularization. From VME04, we know 
that models without regularization give the 
full range of parameters allowed by the data. 
Regularization mainly helps to reduce the sen- 
sitivity of the solutions to noise in the data by 
requiring a smoother sampling of the orbit li- 
braries used in the solution. However, when 
the unregularized solutions give a poor fit to 
the data, regularization does not improve the 
constraints obtained on the parameters. 

Without regularization, the number of or- 
bits used in the fits is generally equal to the 
number of constraints. In the edge-on model, 
orbit libraries employing half the number of 
orbits produced no difference in the results. In 
contrast (as described below), small changes 
in the orbit library size for the inclined model 
indicated some sensitivity to the number of 
orbits used, although not in a systematic man- 
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ner. Since the number of orbits available 
influences the uncertainties in the solution 
(VME04), we analyzed both the inclined and 
edge-on cases with 8100 orbits, providing the 
most uniform comparison of the models. 

4. Results of 3-Integral Modeling 

We present the results of 3-integral model- 
ing with orbit libraries consisting of 8100 or- 
bits each, which were fit to all of the avail- 
able kinematic and photometric data. Fig- 
ure [8] shows 2-dimensional contour plots of 
the total as a function of Mbh (abscissa) 
and Tr for models which assumed that the 
bulge is viewed edge-on (i = 90°). The 
models covered 12 values of Mbh and 8 val- 
ues of Tr (indicated by the grid of points); 
the inner 5 contours are plotted at A(x^) = 
2.30, 4.61, 6.17, 9.21, 11.8 which correspond to 
68.3%, 90%, 95.4%, 99%, 99.73% confidence 
intervals, assuming two degrees of freedom 
(DoF). The innermost contour gives a la up- 
per limit of 4.2 x lO'^ Mq (for 2 DoF) but 
no best-fit value of the BH mass, since the 
values of are all comparable at smaller 
masses. The value of the mass-to-light ratio 
corresponding to the mass upper limit is Tr 
= 2.3, roughly consistent with the estimate 
from the ground-based photometry in § 12.2.21 
Figure [9] shows a 1-dimensional cut through 
the 2-dimensional contour plot at the value 
of Tr = 2.3. The horizontal dot-dash line is 
at Ax^ = 2.30, which corresponds to the la 
confidence interval for two DoF. The vertical 
dashed line in Figure [9] is at the location of 
the best-fit BH mass obtained from reverber- 



ation mapping stud ies: 4:.57~^q^ X 10' Me 
(jBentz et al.ll2006bl '). 

Figure shows a similar 1-dimensional 
cut (at the best-fit Tr = 3.5) through the 
corresponding 2-dimensional space for the 
models that assumed the bulge is viewed at 
the same inclination as the large-scale galaxy 
disk (as inferred from the shape of the disk, 







i = 23°). As an additional test of numerical 
effects, we ran several versions of the inclined 
model with orbit libraries of slightly different 
sizes (ranging between 7600 and 8100 orbits), 
for which we also show the 1-dimensional cuts 
in Figure [TOj There is clearly a minimum in x^ 
between 4 x 10'' and 5 x 10^ M0, but the pre- 
cise location is uncertain. Furthermore, the 
non-smooth nature of the x^ topology sug- 
gests that it is unwise to assign much meaning 
to the formal la x^ intervals since the location 
of the minimum varies quite strongly with or- 
bit library. The most conservative statement 
that can be made is that the models give a 
best-fit BH mass between 4 — 5 x 10^ Mq. 

However, there are two reasons why we are 
reluctant to accept this as more than prelim- 
inary evidence of a "BH detection" based on 
stellar dynamical modeling. First, the best-fit 
X^ value for the inclined model is a factor of 3 
larger than the lowest value of x^ for the edge- 
on model (with the upper limit). Formally, 
therefore, the solution with an edge-on bulge 
and an upper-limit Mbh of 4.2 x 10^ Mq is 
preferred. Second, there is considerably more 
noise in the topology of the 2-D x^ surface 
for the inclined models than for the edge-on 
models, which is indicative of the fact that 
the inclined models are significantly more in- 
fiuenced by the noise in the data than are the 
edge-on models. 

Figure [TT] shows fits for 4 different models 
to the kinematical data from the two ground- 
based datasets (shown by open circles with 
la error bars). The black lines are the fits 
obtained for an edge-on bulge (the solid black 
line is for the model at the upper limit, and 
the dashed black line is for an edge-on bulge 
with no black hole). The red lines are for 
the models with inclination angle i = 23°. 
The solid red line is for the model that gives 
the best fit to the kinematical data, while the 
dashed red line is for the inclined bulge with 
no BH. 
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It is clear that in the case of the edge- 
on models (black lines), the upper limit and 
the "no BH" models provide essentially iden- 
tical fits to the data. For an inclined bulge, 
the model with no BH systematically under- 
predicts the velocity (V), which is why the 
larger BH mass is needed to fit the kinemati- 
cal data. 

If we make the plausible assumption that 
the large-scale inclination of the disk and its 
symmetry axis are likely to be shared by the 
bulge, one can reasonably conclude that the 
best-fit BH mass is 4 — 5 x 10^ Mq, in ex- 
tremely good agreement with the estimate ob- 
tained from reverberation mapping. Models 
in which the bulge grows via secular processes 
or by dissipative minor mergers of disk galax- 
ies p redict such alignment of symmetry axes 



(e.g., Robertson et al. 20061 ) 



We propose this conclusion with some cau- 
tion, however, because the ground-based kine- 
matical data are of relatively low quality and 
do not resolve the sphe re of influence of the 



have shown that it is essential to have large- 
scale two-dimensional kinematical data from 
integral field spectroscopy in order to obtain 
reliable constraints on the inclination angle of 
an ellipsoidal galaxy model from the kinemat- 
ical data alone. 

The existing data and the dynamical mod- 
els presented above show that there is tenta- 
tive, but good agreement between the mass of 
the BH inferred from stellar dynamical mod- 
els and the mass obtained from reverberation 
mapping. 

5. Discussion 

We have presented one of the first direct 
stellar dynamical BH mass estimates for a 
broad-line AGN. We found that if the bulge 
of NGC 4151 was assumed to be viewed edge- 
on, the data only gave an upper limit on the 



dynamical BH mass of 4.2 X 10^ Mq. Al- 
ternatively, when the bulge was assumed to 
have the same inclination to the line of sight 
as the rotation axis of the large-scale disk, 
we obtained a best-fit value for the mass of 
the central BH of Mbh = 4 - 5 x 10^ M©. 
The only other stellar dynamical measure- 
ment of a broad-line A GN is from a recen t 
study of NGC 3227 by IPavies et al.l jioO^), 



who estimated a la range for the BH mass of 
7 - 20 X 10^ M0. The value of Mbh derived 
from stellar kinematical data and 3-integral 
modeling is remarkably consistent with the 
value derived from reverberation mapping in 
both cases. 



BH. Previous studies (jVerolme et al.l l2002l ) |Onken et al 



In addition, iHicks &: MalkanI (j2007l ) have 
analyzed the dynamics of the molecular hy- 
drogen gas in the central regions of both 
NGC 3227 and NGC 4151, and derive BH 
mass estimates for these galaxies that agree 
well with the reverberation mapping esti- 
mates. All of this supports the empirical cal- 
ibration of re verbe ration-based masses from 
jiooi), and indicates that the 
local AGNs used to establish that calibration 
likely already lie very close to the Mbr—ct^ re- 
lation defined by quiescent galaxies. However, 
there is a fair amount of scatter in the AGN 
Mbh — relation, and the mass predicted by 
the relation for NGC 4151 is « 5 x 10^ M© 
(using cj* = 98 km s~^ and the fit to the 



relation by Ferrarese Sz Ford 20051 ). roughly 
an order of magnitude lower than suggested 
by either the reverberation mapping studies 
or the stellar dynamical data presented here. 
While one could attribute the discrepancy of 
the reverberation-based mass in NGC 4151 
to object-to-object variations in the rever- 
beration scaling factor /, a higher value of 
Mbh is supported by the stellar dynamical 
estimate, the gas dynamical estimate, and 
the relationship between BH mass, bolomet- 
ric luminosity, and X- ray variability timescale 
(|McHardv et al.lboOfil ). We note that our stel- 
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lar dynamical models for NGC 4151 are based 
on kinematical data that do not fully resolve 
the sphere of influence of the BH and that 
our modeling is simplistic in that it assumes 
an axisymmetric bulge and takes no account 
of the non-axisymmetry arising due to the 
presence of the galaxy's bar. Thus we con- 
sider this measurement of the BH mass to be 
preliminary. 

We believe that NGC 4151 is an excellent 
target for future stellar dynamical measure- 
ments with higher spatial resolution. One im- 
provement over our present approach would 
be to expand from two slit positions to many, 
or to use an integral field unit (IFU). In the 
seeing-limited regime, observations with the 
optical IFU, SAURON, have re cently been 
published by iDumas et al.l (120071 ) and show a 
similar decrease in the stellar velocity disper- 
sion of NGC 4151 at small radii. While the 
current shutdown of STIS precludes deeper 
studies of the CaT with HST, observations 
of the CO bandhead stellar absorption fea- 
tures at 2.29 ^ m (as was used in the study 
of NGC 3227 bv lPavies et~IDl200fil ) opens the 
possibility of ground-based observations with 
adaptive optics (AO) on large-aperture tele- 
scopes. It is worrying that observations of 
NGC 4151 taken during perio ds of both high 
AGN flux dlvanov et al.ll200d l and low AGN 
flux (|Hicks k Malkanll2007l ) fail to clearly re- 
veal CO bandhead absorption from the cen- 
tral regions of the galaxy, but deeper observa- 
tions may improve the situation. Incorporat- 
ing these types of data into future modeling 
will be a significant step forward. 

The use of high-spatial-resolution AO ob- 
servations with infrared IFUs in conjunction 
with non-axisymmetric modeling will allow us 
to constrain the complex velocity fields aris- 
ing from the bar; this is likely to place even 
stronger constraints on the mass of the cen- 
tral BH. Ultimately, the insensitivity of stellar 
orbits to the outflows and bulk motions that 



affect gas dynamics means that stellar dynam- 
ical BH masses will likely act as the linchpin 
in calibrating the BH masses produced by re- 
verberation mapping. 
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Fig. 1.— 20.5" X 24.5" portion of the HST 
ACS/HRC F550M image of NGC 4151. The 
logarithmic contours atop the Unear grayscale 
show both the bright nuclear point source and 
the lower contrast features at larger radius. 
The indicated scale denotes 1 arcsec. The 
HRC occulting finger can be seen at the top 
of the frame. 




Fig. 3.— Top: {V - R) color in each MDM 
image pixel function of NGC 4151 galac- 
tocentric radius (in arcseconds). Bottom: 
{B — R) pixel color as a function of radius. 




Fig. 2.-— 4.5' X 5.5' portion of the MDM R- 
band image of NGC 4151. The logarithmic 
contours trace the bulge and the nuclear point 
source, while the linear grayscale enhances the 
surface brightness features of the bar. The 
ACS/HRC field of view is overlaid as the thick 
box. The indicated scale denotes 30 arcsec. 



Fig. 4. — Logarithmic contours of the AGN- 
free MGE model of the surface brightness. 
Left: On the size scale of the HST image (as 
in Figure [T|). Right: On the size scale of the 
MDM image (as in Figure [2|). 
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Fig. 5. — i?-band surface brightness along the 
major axis. The AGN-free profile determined 
by the MGE routine and used in the dynam- 
ical modeling is shown as the solid line. The 
raw image profile (including the AGN) and 
the PSF are shown for the HST data (filled 
points and dotted line) and for the MDM data 
(open points and dashed line). 
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Fig. 6. — Kinematic data from KPNO {lef£) 
and MMT {right) parameterized as Gauss- 
Hermite coefficients V , a^hz, and /t4, and the 
rms velocity {-s/V^ + cr^) as functions of spec- 
troscopic slit position. 
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Fig. 7. — Historical light curve of optical con- 
tinuum flux in NGC 4151. The dates of an 
earlier STIS program that used the occulting 
bar (GO-7350), our ToO observations (GO- 
9849), and when STIS stopped functioning are 
indicated. 
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Fig. 8. — Contours of constant derived by 
fitting the kinematical data with 3-integral 
axisymmetric models for a bulge that is as- 
sumed to be viewed edge-on. The dots corre- 
spond to the values of parameters (Mbh^T^r) 
for which solutions were obtained. The three 
inner most contours (moving left to right at 
R~ 3) correspond to the la, 2o", and 2>a con- 
fidence intervals respectively (68.5%, 90% and 
95%). The contour plot is obtained using a 
smoothing kernel. 
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Fig. 9. — values as a function of Mbh ob- 
tained by taking a cut through Figure [8] at 
Tr = 2.3. The horizontal dot-dash line is 
at Ax^ = 2.3 (68.5% confidence interval for 
2 DoF) above the minimum value of x^- The 
vertical dashed line indicates the best-fit value 
obtained from reverberation mapping. 
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Fig. 10. — x^ values as a function of Mbh 
obtained by taking cuts at Tr =3.75 through 
several dynamical models with inclined bulges 
that are created with slight variations in the 
number of kinematic orbits (7600-8100). The 
vertical dashed line indicates the best-fit value 
obtained from reverberation mapping. 
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Fig. 11. — Fit to kinematical data 
parametrized as Gauss-Hermite coefficients 
V, a, /13 , /i4 as a function of position along the 
slit |x| for ground-based data from KPNO 
(top four panels) and MMT (bottom four pan- 
els). Solutions used libraries of 8100 orbits 
for 4 different values of Mbh as indicated by 
the labels. The black curves are fits obtained 
with two edge-on models: Mbh = 4 x 10'' Mq 
(solid line) and Mbh = (dashed line). The 
red curves are for models with inclination an- 
gle i = 23°: Mbh = 4.5 x lO'^ Mq (solid line) 
and Mbh = (dashed line). 
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Table 1 

AGN-SUBTRACTED MuLTI-GaUSSIAN EXPANSION RESULTS 



Gaussian 


Surface Density 


Gaussian Sigma 


Axial 


Number 


(L0,i? 


(arcsec) 


Ratio 


1 


8.75E+04 


0.14 


1.00 


2 


2.43E+04 


0.32 


1.00 


3 


8.38E+03 


1.19 


1.00 


4 


3.19E+03 


2.44 


1.00 


5 


1.25E+03 


6.46 


1.00 


6 


2.75E+02 


11.57 


1.00 


7 


1.50E+02 


45.21 


0.67 


8 


5.06E+00 


203.53 


1.00 


9 


0.88E+00 


203.53 


0.27 



Table 2 

Gauss-Hermite Fits to Kinematic Data 



Position 



Dataset 


Along Slit^ 
(arcsec) 


Velocity 
(km s~^) 


a 

(km s"-*^) 




hi 


MMT 


12.0±0.3 


-10±11 


76±13 


O.OlibO.ll 


0.00±0.12 




11.4±0.3 


-14±8 


46±11 


0.02±0.14 


0.00±0.18 




10.8±0.3 


-18±8 


37±11 


-0.02±0.14 


-0.02±0.18 




10.2±0.3 


-11±8 


45±10 


0.00±0.13 


0.03±0.17 




9.6±0.3 


-26±8 


63±10 


O.OOibO.lO 


0.05±0.11 



^Uncertainties indicate aperture size along the slit. 

Note. — The kinematic data were symmetrized in the subsequent dynam- 
ical modeling. [The complete version of this tabic is in the electronic edition 
of the Journal. The printed edition contains only a sample.] 
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